Shock waves and population inversion in collisions of ultracold atomic clouds 
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Ultracold atomic gases represent an ideal toolbox to study quantum effects that are difficult to 
probe using other systems. Here, we use a time-dependent density matrix renormalization group 
approach to show that the collision of two interacting bosonic clouds in one dimension gives rise 
to shock waves with a concomitant local energy distribution typical of population inversion, i.e., 
an effective negative temperature. A classical hydrodynamic description compares well with the 
exact quantum dynamics only up to the gradient catastrophe time. Such a highly nonequilibrium 
local distribution, however, does not prevent the system from recovering its initial state after an 
oscillation period which is renormalized by the interaction. All these results can be easily tested 
experimentally. 
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Experiments on ultracold atomic gases have unveiled 
new aspects of quantum systems comprising a large num- 
ber of particles, i.e., of quantum fields. The degree of iso- 
lation from the external environment and the long time 
scales allow to study details of the quantum dynamics 
that are not easily accessible, e.g., in solid-state systems. 
Moreover, the dynamics are often nonlinear and a sub- 
stantial effort is needed to restrict the system to a regime 
where linear response theory is applicable [1]. As a con- 
sequence, the concept of quasiparticle - an elementary 
excitation above the ground state - is hardly of any help 
in understanding the time evolution. Instead, at very low 
temperatures the system behaves rather as a single col- 
lective object, and a description in terms of fields rather 
then point particles is more appropriate. 

Examples of collective quantum dynamics are collision 
experiments between clouds of ultracold atomic gases. 
Various phenomena have been studied in this context 
such as the interference of matter waves [2], dispersive 
shock- waves in BECs [3, 4], superfluidity, shock- wave for- 
mation and domain wall propagation in the unitary Fermi 
gas [5, 6], spin transport [7], and lack of thermalization 
in quasi-integrable ID systems [8, 9]. Most recently [10] 
a negative temperature, namely a population inversion in 
the energy distribution of the motional degrees of free- 
dom of atomic gases at equilibrium has been realized. In 
this work it was also noted that such a population inver- 
sion does not necessarily imply a fast decay to the true 
thermal equilibrium state, thus showing the quite unique 
properties these systems possess. 

In this Letter we show that out of equilibrium distri- 
butions characterized by an effective negative temper- 
ature have also been realized in the "quantum New- 
ton cradle" of Ref. [8] for almost any value of the in- 
teraction strength. With reference to this experimental 
set up, we study the collision of two clouds of one di- 
mensional bosons for arbitrary interaction strength, by 
means of a time-dependent density matrix renormaliza- 
tion group (TDMRG) approach [11-14]. We demonstrate 
both the emergence of population inversion (negative ef- 



fective temperature) at the onset of shock waves - where a 
hydrodynamic description deviates from the exact quan- 
tum dynamics - and no visible decay in the density profile 
or in the local energy distribution function as a func- 
tion of time, despite the presence of population inver- 
sion. Instead, the dynamics are essentially periodic with 
the oscillation frequency renormalized by the interaction 
strength. All these results can be easily verified experi- 
mentally. 

Model and methods — The Lieb-Liniger model [15, 16] 
provides an excellent description of one-dimensional ul- 
tracold bosonic atoms [17, 18]. In terms of a bosonic field 
^/(x) its Hamiltonian reads 
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where gs G [0, +oo] is a coupling constant and m is the 
atom mass. In the following we will consider a time- 
dependent external potential V(x,t) changing abruptly 
at t = (quench protocol). Hamiltonian (1) is integrable 
for any gg when V(x) = 0, while for V(x) ^ the ex- 
act eigenstates and eigenvalues are known for free bosons 
gs = and hardcore bosons gg = +oo, the latter being 
equivalent to free fermions (Bose-Fermi mapping theo- 
rem [19, 20]). 

In the following we consider a closed finite system of 
N = 20 particles [21] initially prepared in the ground 
state of a double well potential of the form 
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At subsequent times t > the trapping potential be- 
comes 



V(x,t > 0) = \rLu 2 x 2 ■ 



(3) 



and the two clouds are forced to collide (microcanonical 
picture of transport [22, 23]). 
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FIG. 1. (Color online) a) Snapshots of the density profile p(x,t) in units of 1/a at various times for gs = Q.2J/a (7 = 2). 
The black lines are the exact TDRMG results, the grey lines are obtained from the solution of Eq. (4). At t — the two 
profiles are indistinguishable, but they deviate from each other at the onset of shock-wave formation (t = 0.38r). b) Energy 
distribution function f(x, e, t) (6) calculated from the exact Wigner function (5) obtained with TDMRG. The color plots show 
the values of f(x,s,t) in the (x,e) plane at the corresponding times in panel a). The vertical black, grey and light grey lines 
correspond to the sections for fixed x of f(x, e, t) shown in Fig. 2. Note that shock-wave formation is characterized by a highly 
nonequilibrium distribution. c)-d) Same as panels a) and b) but for gs = 2.0J/a (7 = 20). The shock-wave for gs = 2.0J/a 
forms at roughly the same time t = 0.32r as for gs = Q.2J/a. 



We are able to access the exact dynamics of (1), (2) 
and (3) with a TDRMG approach. TDRMG has been ap- 
plied mainly to lattice systems for relatively short time 
scales, but simulations of systems in the continuum limit 
and for quite long time scales (of the order of several 
trapping potential periods) are feasible [26-30]. In our 
simulations we have employed two different discretiza- 
tions of (1) [26, 31], either as a Bose-Hubbard model 
(nonintegrable), or as a XXZ spin chain (integrable) us- 
ing the Bose-Fermi mapping for arbitrary gs [24-26]. No 
substantial difference in the results has been observed. 

In the following, lengths are expressed in units of the 
lattice spacing a of the discretized model, which is a 
small, but otherwise arbitrary length scale, energy is 
in units of J — h 2 /(2ma 2 ), and time in units of the 
post-quench oscillation period r = 7r/w2- The interac- 
tion parameter gs is given in units of J/a. Occasionally 
we use the Lieb-Liniger parameter 7 = mgs / {ti 2 p) with 
p = 0.05/a, an indicative value of the density in the inho- 
mogeneous system considered here. In our simulation we 
used a lattice of length L = 600a and D — 120a [see (2)]. 
We ensured, by separately tuning uji and W2 for each 
value of the coupling gs, that the particle density per 
site never exceeds 0.15, thus lattice effects are negligible 
(continuum limit) [28, 29, 31]. 

Numerical results and hydrodynamic description - 
TDMRG is based on the matrix product state repre- 



sentation of the full wave function [14]. However the 
Runge-Gross theorem [32] of Time-Dependent Density 
Functional Theory (TD-DFT) guarantees that an ex- 
act hydrodynamic description (purely in terms of density 
and velocity fields) of quantum dynamics exists [22] , al- 
though the analytical expression of the stress tensor is 
unknown even for free fermions [33, 34]. To our knowl- 
edge the best possible hydrodynamic description for the 
present case is a generalized nonlinear Schrodinger equa- 
tion (GNLSE) [18, 35-38] 

h 2 

ihd t ^(x,t) = -—8%V(x,t)+(l>(j>)V(x,t)+V(x,tW(x,t). 
2m 

(4) 

where &(x) is a complex field, V(x,t) is specified by (2) 
and (3), p(x) = ^(x)] 2 is the density, and <fi(p) the Gibbs 
free energy per particle which can be obtained in the lo- 
cal density approximation (LDA) from the Bethe Ansatz 
solution of the Lieb-Liniger model. Accurate numerical 
values of this quantity are available [18]. We employ this 
hydrodynamic description to better identify the time of 
formation of shock waves. 

TDRMG and hydrodynamics are compared in Fig. 1 
for <7b = 0.2 J/a and gs — 2.0 J/a. The excellent agree- 
ment for t = indicates that lattice effects are negligi- 
ble (continuum limit) and that LDA works well for just 
N = 20 particles. The hydrodynamic description is seen 
to model well the dynamics of the density immediately 
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FIG. 2. Energy distribution function f(x,e,t) corresponding 
to first (t — Or) and third (t ~ 0.3r) snapshots in Fig. 1, pan- 
els b) and d) . The black, grey and light grey lines correspond 
to x/a = 25, 40, 60 for the two upper quadrants (gs = 0.2J/a) 
[Fig. 1 panel b)], and x/a = 15, 30, 50 for the two lower quad- 
rants (gB = 2.0 J /a) [Fig. 1 panel d)], respectively. Note that, 
neglecting oscillations for large e - due to the finite number 
of particles - the initial distribution decreases monotonically, 
while a maximum develops for finite e after the shock-wave 
formation t ~ 0.3r, i.e. a population inversion. 



after the quench. However it fails when oscillations form 
in the GNLSE solution (t = 0.38r Fig. la, t = 0.32r 
Fig. lc, grey lines), corresponding to the formation of 
shock waves (gradient catastrophe [39, 40]). These shock 
waves with oscillatory behaviour are known as dispersive 
and occur in inviscid fluids, e.g. Bose-Einstein conden- 
sates [3, 4, 40]. 

Our TDMRG results are very similar to the exper- 
imental data reported in Ref. 5 where a viscosity was 
introduced in the hydrodynamic equations to describe 
shock waves, while in the TD-DFT calculation in Ref. 41 
a renormalized kinetic term A<9 2, J' was used for the same 
reason. This is not justified here since Eq. (4) has no free 
parameters and it is an excellent approximation up to the 
gradient catastrophe for any g B . Introducing viscosity 
would contradict the fact that almost no dissipation is 
present in our system (see below). 

Population inversion and negative temperature — In 
order to study in more detail the particle distribution 
during the quantum shock wave, we use the Wigner func- 
tion [33, 34] 



W(x,p,t) = 



1 



dy p(x + y,x - y;t)t 



(5) 



where p(x\x;t) = (W (x',t) \& (x,t)) is the one-body 
density matrix. The exact one-body density matrix is 
easily accessible within TDRMG. Neglecting negative 
values, W(x,p, t) can be thought of as a local momen- 
tum (p) distribution as in the Boltzmann equation [43]. 
The local energy distribution f(x, e, t) is defined with 
respect to the local comoving (Lagrangian) reference 
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FIG. 3. Inverse information compressibility as a function of 
time for various interaction strengths gsa/J = 0.1, 0.5, 1.0. 
Note the divergence in correspondence to a fully developed 
shock wave at t = 0.4-^0.5r. See Fig. 1, snapshots at t = 0.38r 
in panel a) and at time t — 0.32r in panel c). 



frame with velocity mv(x,t) — j(x,t) / p(x 1 t) [22], with 
j(x) = J dppW(x,p,t) and p(x 1 t) = J dpW(x 1 p,t). 
Thus J dpW(x,p — mv(x),t) = and the energy dis- 
tribution reads [42] 

f(x, e, t) — 2irh ^ W(x, s\/2me — mv(x)) . (6) 



This is the quantity shown in the color plots in Fig. 1 and 
for selected values of x in Fig. 2. At t — the distribution 
f(x,e,t) decreases monotonically with e (leaving aside 
oscillations related to the finite particle number) indi- 
cating an equilibrium energy distribution. In correspon- 
dence of the shock-wave formation t ~ 0.3t the energy 
distribution is no longer an equilibrium one, f(x,s,t) is 
larger for values of e away from zero (see Fig. 1 and quad- 
rants at the right of Fig. 2), signalling a population in- 
version, namely an effective negative temperature. Pop- 
ulation inversion leads to the break down of the LDA, 
namely to the departure from the GNLSE (Eq. 4) solu- 
tion. Notice also that we have observed nonequilibrium 
distributions during the evolution for gg > 0.02J/a, a 
rather small value of the interaction parameter corre- 
sponding to an average Lieb-Liniger parameter of 7 = 
0.2, well within the weakly interacting regime, thus sug- 
gesting that such phenomena should be present at all 
interaction strengths. 

In order to corroborate the presence of a negative tem- 
perature out of equilibrium we characterize the state 
of the system with the information compressibility [44], 
which measures the relative change of number of avail- 
able microstates of an open system in response to an en- 
ergy variation. Given a system with density matrix p(t) 
and Hamiltonian H, the energy is E(t) = Tr[p(t)'H] and 
the thermodynamic entropy S(t) = — fc B Tr [p(t) \np(t)]. 
Call f2 the number of microstates available to the sys- 
tem. The information compressibility is defined as the 
relative variation of the number of microstates with re- 
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FIG. 4. (Color online) Upper plot, deviation A(t) as a func- 
tion of time for gs = 0.5J/a. The vertical notches indicate 
the instants (t — 1.09r and 2.18r) where the density is clos- 
est to the initial one. The corresponding density profiles are 
shown in the middle panel (black lines) along with the hydro- 
dynamic results which do not show periodicity (grey). The 
density profiles at different times have been shifted in the ver- 
tical direction. In the lower plot the renormalized oscillation 
period t*(qb)/t is shown, with r = n/uj2, as a function of the 
interaction strength gs extracted from the minima of A(t) as 
in the upper plot. The red triangles show the results ob- 
tained with the XXZ spin chain discretization while the blue 
circles with the numerically more demanding Bose-Hubbard 
discretization. The red line at the bottom left represents the 
frequency extracted from the gs = exact data. It deviates 
from the exact result t* (0) = r since lattice effects distort the 
density profile in time [31]. This is not an error due to the 
approximation introduced by TDMRG. 



spect to the energy variation at time t [44] 
K T (t) 



i sn 

TlSE 



(7) 



Using the microcanonical relation Q — exp (S/k B ) one 
arrives at the computationally more convenient definition 



Kj{t) = 



(8) 



Note the similarity of this quantity with the thermody- 
namic definition of inverse temperature [44]. Here, we 
take for p(t) the density matrix of half of the system, 
i.e., relative to the interval [— L/2, 0]. The results for the 
inverse information compressibility are shown in Fig. 4 
for various interaction strengths. The interesting point 
is the divergence of KJ 1 ^) in all cases for t/r ~ 0.4, i.e., 
in correspondence to a fully-developed shock wave, and 



negative values of this quantity at later times. If we in- 
terpret the inverse compressibility as an effective temper- 
ature, this behaviour is clearly suggestive of a population 
inversion, in agreement with our previous results. 

Oscillation frequency shift — Finally, we note that al- 
though the system explores strongly nonequilibrium en- 
ergy and momentum distributions we find that the initial 
density profile and the initial Wigner function are recov- 
ered after a time of order r for any gs £ [0, +oo]. This 
periodicity is clearly seen in the mean square deviation 
of the density profile at time t from the initial one 



A(*) = 




dx (n(x, t) — n(x, 0)) 2 



(9) 



shown for g B — 0.5 J/a in the upper panel of Fig. 4. 
A(t) essentially drops to zero at times t — 1.09r and 
t = 2.18t indicating that the system has returned to the 
initial state, at least within our simulation times. The 
density profiles corresponding to the minima of A(t) are 
shown in the middle panel of Fig. 4 alongside with the 
results of classical hydrodynamics. While the exact dy- 
namics are essentially periodic, the classical one are not. 
We remind the reader that the static density (as shown 
in the lower t — Or profile in Fig. 4) and the dynam- 
ics up to the gradient catastrophe are well captured by 
the GNLSE (4). A similar analysis can be performed for 
all gB and the renormalized period T*(gs) is shown in 
the lower panel of Fig. 4. In the exactly solvable limits 
gB = and gB = +oo the period is not renormalized. In 
between these extrema it has a nonmonotonic behaviour 
with a maximum in the interval 0.01 < gB a /J < 0.1 
(0.2 < 7 < 2). The almost perfect periodicity observed 
for any gs is a strong indication of very small dissipation, 
in agreement with experimental results [8] . 

Conclusions — Using TDMRG we have computed the 
exact dynamics of interacting bosons in a ID collision 
experiment for arbitrary interaction strengths, and com- 
pared the results with a parameter-free hydrodynamics 
deduced from the exact Bethe-Ansatz solution in the 
LDA. Even for small coupling, shock-wave formation 
leads to the break down of the hydrodynamic descrip- 
tion. We have shown that quantum shock waves lead to 
a population inversion in the local energy distribution, 
namely to a negative effective temperature. Finally, we 
have computed the shift in the oscillation period due to 
interactions, a prediction which can be easily tested in 
experiments. Our results suggest that statistical ensem- 
bles with negative temperatures for the motional degrees 
of freedom, as shown in Ref. [10], are a common feature 
in collision experiments with ultracold gases [8] . 
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